
// @BeginInterface

// Texture access using a bicubic filtering.
//
// @param tex		Texture to be sampled.
// 
// Important: 'tex' magnification/minification filter shall be set to GL_NEAREST.
vec4 textureBicubicFilter(sampler2DRect tex, vec2 coord);

// Texture access using a bicubic filtering (faster).
//
// @param tex		Texture to be sampled.
// 
// Important: 'tex' magnification/minification filter shall be set to GL_LINEAR.
vec4 textureFastBicubicFilter(sampler2DRect tex, vec2 texCoord);

// @EndInterface

vec4 textureCubicFilter(float a, vec4 c00, vec4 c10, vec4 c01, vec4 c11);

vec4 textureBicubicFilter(sampler2DRect tex, vec2 texCoord)
{
	// Center textel coordinate to avoid round-off errors
	texCoord -= vec2(0.5, 0.5);

	float px = floor(texCoord.x);
	float py = floor(texCoord.y);
	float fx = texCoord.x - px;
	float fy = texCoord.y - py;

	return (
		textureCubicFilter(
			fy,
			textureCubicFilter(fx, texture(tex, vec2(px - 1, py - 1)), texture(tex, vec2(px + 0, py - 1)), texture(tex, vec2(px + 1, texCoord.y - 1)), texture(tex, vec2(px + 2, py - 1))),
			textureCubicFilter(fx, texture(tex, vec2(px - 1, py + 0)), texture(tex, vec2(px + 0, py + 0)), texture(tex, vec2(px + 1, texCoord.y + 0)), texture(tex, vec2(px + 2, py + 0))),
			textureCubicFilter(fx, texture(tex, vec2(px - 1, py + 1)), texture(tex, vec2(px + 0, py + 1)), texture(tex, vec2(px + 1, texCoord.y + 1)), texture(tex, vec2(px + 2, py + 1))),
			textureCubicFilter(fx, texture(tex, vec2(px - 1, py + 2)), texture(tex, vec2(px + 0, py + 2)), texture(tex, vec2(px + 1, texCoord.y + 2)), texture(tex, vec2(px + 2, py + 2)))
		)
	);
}

float w0(float a)
{
//	return (1.0f/6.0f)*(-a*a*a + 3.0f*a*a - 3.0f*a + 1.0f);
	return (1.0f/6.0f)*(a*(a*(-a + 3.0f) - 3.0f) + 1.0f);
}

float w1(float a)
{
//	return (1.0f/6.0f)*(3.0f*a*a*a - 6.0f*a*a + 4.0f);
	return (1.0f/6.0f)*(a*a*(3.0f*a - 6.0f) + 4.0f);
}

float w2(float a)
{
//	return (1.0f/6.0f)*(-3.0f*a*a*a + 3.0f*a*a + 3.0f*a + 1.0f);
	return (1.0f/6.0f)*(a*(a*(-3.0f*a + 3.0f) + 3.0f) + 1.0f);
}

float w3(float a)
{
	return (1.0f/6.0f)*(a*a*a);
}

vec4 textureCubicFilter(float a, vec4 c00, vec4 c10, vec4 c01, vec4 c11)
{
	return (c00 * w0(a) + c10 * w1(a) + c01 * w2(a) + c11 * w3(a));
}

// g0 and g1 are the two amplitude functions
float g0(float a)
{
	return (w0(a) + w1(a));
}

float g1(float a)
{
	return (w2(a) + w3(a));
}

// h0 and h1 are the two offset functions
float h0(float a)
{
// note +0.5 offset to compensate for CUDA linear filtering convention
	return (-1.0f + w1(a) / (w0(a) + w1(a)) + 0.5f);
}

float h1(float a)
{
	return (1.0f + w3(a) / (w2(a) + w3(a)) + 0.5f);
}

vec4 textureFastBicubicFilter(sampler2DRect tex, vec2 texCoord)
{
	// Center textel coordinate to avoid round-off errors
	texCoord -= vec2(0.5, 0.5);

	float px = floor(texCoord.x);
	float py = floor(texCoord.y);
	float fx = texCoord.x - px;
	float fy = texCoord.y - py;

	float g0x = g0(fx);
	float g1x = g1(fx);
	float h0x = h0(fx);
	float h1x = h1(fx);
	float h0y = h0(fy);
	float h1y = h1(fy);

	return (
		g0(fy) *
		g0x * texture(tex, vec2(px + h0x, py + h0y)) +
		g1x * texture(tex, vec2(px + h1x, py + h0y)) +
		g1(fy) *
		g0x * texture(tex, vec2(px + h0x, py + h1y)) +
		g1x * texture(tex, vec2(px + h1x, py + h1y))
	);
}